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ABSTRACT 

The equilibrium and stability properties of a coupled two-component BEC is studied us- 
ing a variational method and the one-dimensional model of Williams and collaborators. The 
variational parameters are the population fraction, translation and scaling transformation of 
the condensate densities, assumed to have a Gaussian shape. We study the equilibrium and 
stability properties as a function of the strength of the laser field and the traps displacement. 
We find many branches of equilibrium configurations, with a host of critical points. In all 
cases, the signature of the onset of criticality is the collapse of a normal mode which is a 
linear combination of the out of phase translation and an in phase breathing oscillation of 
the condensate densities. Our calculations also indicate that we have symmetry breaking 
effects when the traps are not displaced. 
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I. INTRODUCTION 



One of the most interesting achievements in the field of boson condensation was the 
experimental observation of binary mixture of trapped condensates [1-4]. The species in the 
mixture can be atoms with different F spin orientations [3-4] or simply different hyperfine 
states of the same atom [1-2]. In the last case we can have inter-conversion between the 
components through the coupling of the atoms with an external laser field [5]. 

In the case of uncoupled components many effects have been predicted theoretically 
and determined experimentally such as spatial phase separation [6], stability properties of 
the equilibrium states [6-7] and symmetry breaking effects [6,8]. These aspects of binary 
condensate mixtures have been treated in the literature using various methods such as the 
Thomas- Fermi (TF) approximation [8-9] and numerical solutions of the appropriate coupled 
Gross-Pitaevskii (GP) equations [10]. 

In this paper we use the variational method [11-13] to study the equilibrium properties of 
a coupled two-component BEC. In our exploratory calculation we use the one-dimensional 
model of reference [14] and take as variational parameters the population fraction, translation 
and scaling transformation of the equilibrium state densities, all assumed to have a Gaussian 
shape. We consider the intensity of the interaction between the condensates equal and the 
detuning is put equal to zero, leading to equal equilibrium population fraction. In our 
calculation we investigate the behavior of the spatial phase separation as a function of the 
trap displacement and the strength of the laser field. The coupling with the laser field have a 
stabilizer effect in the process of phase separation, opposite to the one coming from the trap 
displacement. Our paper is organized as followsdn section 2 we show briefly how to use the 
variational method to study the equilibrium and stability properties of a binary mixture of 
condensates. Specifically, we derive the general form of the equations of motion and develop a 
scheme to find the normal modes. In section 3 we present our numerical results and discuss 
its physical significance. A summary of our numerical results is presented in section 4. 

II. VARIATIONAL STUDY OF THE EQUILIBRIUM AND STABILITY 
PROPERTIES OF A BINARY MIXTURE OF CONDENSATES. 
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A. The variational approach 



The starting point of our discussion is the action [15] 

S = J dzdt[J2^' l Pj(z,t)ip j {z,t)-£(z,t)] (1) 

3 

where ipj(z,t), j — a,b are the condensate wave function of each component in the mixture 
and €(z,t) is the energy density, 

e(z,t) = T,^[-^-f- 2 + vl p ]^ + \ E *kjWm 2 + n[^6 + r b 4>a] . (2) 



3 k,j 

In the above expression Vf rap (z), j — a,b are the trapping potential of each component 

Vlap = + ljZ ) 2 (3) 

where 7 a = — 1 and 7& = 1, zq is the trap displacement, \ aa = and A a & are the strength 
of the intraspecies and interspecies interactions, respectively, and the last term comes from 
the coupling with the laser field responsible by the interspecies tunnelling, with Q the Rabi 
frequency. 

Imposing that the action (1) is stationary with respect to a variation of ipj(z,t) subject 
only to the normalization constrain \^j( z ,t)\ 2 — N> leads to the coupled time-dependent 
Gross-Pitaevskii equations for the condensate wave functions [15], 

= ( _ ^^ + v *& + Xaa ^ 2 + A <* W 2 )^« + ( 4 ) 

in w = + Vt ™ p + xm2 + + ^ (s) 

To get the equations for the stationary states, we look for solutions of equations (4) and (5) 
of the form ipj(z,t) = e~ l ^ipj{z) , which gives rise to the time-independent GP equations: 

h 2 d 2 

M» = (~2^^2 + V trl + X aa\H 2 + A a6 |^| 2 )^a + Wb (6) 

h 2 d 2 

H4>b = (~2^^2 + V *°P + ^ab\^a\ 2 + MM^b + Jtya (7) 

We can view the two condensate wave functions as components of a spinor of "quasi-spin" 
equal to | [16] which leads immediately to the property that the stationary equations are 



invariant under a transformation which is the product of a space reflection, P z , and an atom 
exchange, P eX) where P ex = iexp — %\a x . Therefore, as (P z P ex ) 2 = P z Pex, we can classify the 
stationary states as even (gerade) and odd (umgerade) under this transformation. In the 
even class ip%(z) = ip%(—z), whereas in the odd class VCl- 2 ) = — ^bi~ z )- The next step is to 
solve numerically the coupled stationary equations for the condensate wave functions to map 
the solutions as a function of Q and z . However, before embarking on this task, we found it 
worthwhile to adopt a simpler approach, the variational method. In the variational method 
the search for stationary states reduces to finding the stationary points in a finite dimensional 
energy surface which is much simpler than the corresponding search in an infinite dimensional 
space, when we solve exactly the coupled GP equations. Besides, the variational solutions 
can be used as an initial guess in the numerical solution of the coupled GP equations [10]. 

In the variational approach we parametrize the time dependence of the condensate wave 
function through a set of 2d parameters [12-13] which we denote by w = (wi, w 2 , . . . , w 2 d), 

^(z,t)=^(z,w(t)) (8) 

When we replace the condensate wave functions parametrized as in equation (8), in equa- 
tion (1), the action reduces to a "classical" action, whose variation leads to Hamilton-type 
equations of motion in terms of these parameters [18] 

Er«(w)w.= ^— (w), ( 9 ) 

l OWk 

where E(w) is the spatial integral of the energy density (2), with the condensate wave 
functions parametrized as in equation (8) 

E(w) = [ dz£(z,w) (10) 



and the antisymmetric matrix r^(w) is given by 

r H (w) = -2I m S/^M(,,w)^,w) 

It follows from equation (9) that the equilibrium configurations are determined by the equa- 
tions 

f)P 

^_(w°) = 0, fe = l,2,...,2d. (12) 
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B. Normal Modes 



To investigate the stability of the equilibrium configurations, we calculate the energies of 
the normal modes. They are stable if the energies are real and positive and unstable if one 
of the energies is complex. To find the normal modes, we linearize the equations of motion 
in the neighborhood of an equilibrium configuration, leading to: 



(13) 



In this equation w are the displacements from equilibrium w = w° + w and H^, Tki are, 
respectively, the Hessian and the matrix T(w) evaluated at the equilibrium configuration, 
that is 

d 2 E , n , 



H 



ki 



dw k dw t 



w 



and 



r kl = r fcZ (w°) 



(14) 



(15) 



To proceed, we divide the 

2d parameters into two groups, w = (gi, q 2 , . . . , qd,Pi,P2, ■ ■ ■ ,Pd) = (q, p), of coordinates 
q and momenta p. Schematically, the reasoning behind this splitting is that the amplitude 
of the condensate wave function, whose square is the condensate density, depends only on 
the coordinates while its phase, whose gradient is the velocity field, depends on the momenta 
and the coordinates. 

In terms of the coordinates and the momenta the equations of motion (13), read 



(16) 



where T and H are, respectively, antisymmetric and symmetric matrices which can be written 
in terms of four d x d blocks 

/ 





= H 











( 



\ 



V 



qq 



pq 



qp 



H 



Hqq H. 



\ 



qp 



■ pp 



J 



V 



H pq H 



pp 



(17) 



where, for example, T qq and H qq are d x d matrices whose elements are given by equations 
(11) and (14), where the derivatives are with respect to the coordinates, with an analogous 
definition for the other d x d matrices. 

In our method to find the normal modes, we try to stay as close as possible to the one 
adopted in the standard case [17]. To begin with, we should find a transformation to a set 
of new coordinates and momenta 



(18) 





= T- 1 


M 









such that 



r -i HT i 



O A 

-A O 



(19) 



where A is a diagonal d x d matrix whose diagonal elements are the normal mode energies. 
In terms of the new coordinates and momenta, the equations of motion, (16), reduce to 



A 

-A 



(20) 



leading to 



Qk — A fc P fc , P k — —A k Q k 



(21) 



The equations (18)-(19) define the transformation to the normal modes and our scheme 
to find the matrix T _1 is a straightforward generalization of the standard case [17]. The 
starting point is to solve the eigenvalue problem 



3 -i H y(fc) = AV (k) 



(22) 



where S is the hermitian and antisymmetric matrix S = — iT 

As in the standard case, this eigenvalue problem has two properties: (i) if is an 
eigenvector with eigenvalue A fc , then V^* is also an eigenvector with eigenvalue — A* k . (ii) 
The eigenvectors with different eigenvalues are S orthogonal, that is yW + 3V'( fc ) = 0, if 
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If one of the eigenvalues of equation (22) is complex, the system is unstable. If all the 
eigenvalues are real it is stable and we can find the transformation to the normal modes 
as follows. We define a matrix S -1 , whose first d columns are the 2d components of the d 
eigenvectors with real positive eigenvalues and the next d columns, the 2d components of 
the corresponding eigenvectors with real negative energies. 

In terms of S _1 , the eigenvalue equations (22) reads 



S 'HS 



A 








(23) 



and the orthogonality of the eigenvectors gives that 

/ 



s-^ss- 1 

The matrix T -1 is related to S" 1 by 

T -i 

where U is the unitary matrix 



U = -P 



1 

-1 



(24) 



S _1 U 



(25) 



1 i 

1 -i 



(26) 



To summarize, we state the main steps in our procedure: (i) First we determine the 
equilibrium configurations by solving the 2c? equations (12). (ii) Next, for each equilibrium 
configuration, we solve the eigenvalue equation (22). (iii) When the configuration is stable, 
we construct the matrix S _1 from the eigenvectors as indicated above and T _1 as shown 
in equation (25). When the system oscillates in a normal mode only one pair (Qk,Pk) is 
different from zero and equation (18) gives the time evolution of the system in this case. 

C. Gaussian ansatz 

Now that we have established the general framework of our calculations, we turn to our 
specific application where the variational parameters are related to the population fraction, 
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translation and breathing shape oscillation of the condensate densities. Thus, the condensate 
wave functions are written as [18] 

ipj(z, w) = e iF ^ z ' w) Aj(z, w) (27) 

where the amplitude of the condensate wave function is parametrized as 

1 ~^7~ 



Aj(z,w) = JNnj --e 2q ^ (28) 



1/2 



q2j 



and the phase as 



Fj(z, w) = Plj (z - Qlj ) + £L(z- q l3 f + 9 3 (29) 

Zq 2 j 



The parameters introduced in the above expression have the following physical interpretation: 
rij is the population fraction of each condensate, 

N Uj = J \ipj{z,w)\ 2 dz , (30) 

qij is the center of mass of the spatial distribution of each condensate, 

Nrijqij = J z\^j(z,w)\ 2 dz , (31) 

and q 2 j/V2 is the width of the spatial distribution of each condensate 

2 

N Uj q f = J(z- qij ) 2 \^(z,w)\ 2 dz, (32) 

The momenta pij,p 2 j and 9j appear in the phase of the condensate wave functions. The 
Pij is the center of mass moment of each condensate 

NnjPlj = -iJ ^ w)^(z, w)dz (33) 

and p 2 j is connected to the expectation value of the dilatation operator in the center of mass 
frame of each condensate 

^Nn jP2j q 2j = J (tp*(z, w)^(^, w) - w)^-(z, w)) (z - q Xj )dz , (34) 

showing that qij,Pij is related to the translational degrees of freedom and q 2 j,p 2 j to the 
breathing shape oscillation. 



Since E(w) depends only on the phase difference between the two condensates, the 
equations of motion (9) give that the total number of particles is conserved, n a (t) + rib(t) = 1. 
This reduces the number of degrees of freedom to ten, which we take as the relative phase 9 = 
9b — 9 a and the relative population fraction n = n fr~ n ° , in addition to the q±j, q2j , Pij , P2j , j = 
a, b. 



To perform the calculations we specify the model parameters . Following reference [14], 
we consider our system to be 87 Rb and take JV = 2.3 x 10 4 , v z = 60Hz which gives z sho = 



yfi/muj z = 1.4/im and we put all the interaction strength equal to X aa = Xbb = X a b = 17.5 
/iin sec _1 fr, in order to reproduce qualitatively the spatial distribution of the condensates 
shown in FIG. 2 of reference [14]. 

Our first task is to solve the ten equilibrium equations (12). Six of the equations lead 
immediately to pj- = p°j = 0, j = a, b, equal equilibrium population fraction n° = 0, and 
9° = or 7i. The solution with 9° = belongs to the even class under P z P ex and the ones 
with 9° = 7r, to the odd class. We restrict the calculations to the odd class since, as Q > 0, 
the lowest energy configuration necessarily belongs to this class. Besides, to be an eigenstate 
of P z P ex the equilibrium parameters should obey the conditions = — q® b and q\ a = q^- 
Therefore, to find the equilibrium configurations we calculate the zeros of four functions 
4^-(w°) = with k — 1,2 and j — a.b and the parameters restricted as indicated in the 



above discussion. 

We characterize the equilibrium configurations (eqc) by the localization of the center of 
mass of each component and in Fig.l we have a graph of the relative distance between the 
centers of mass as a function of Q, for appropriately chosen values of zq. As shown in Fig.l, 
for zq = 0.23z s ho we have only one branch of eqc. For small values of Q/huj z the condensates 
are well separated and when the intensity of the laser field increases the overlap between the 
condensates increases very slowly up to a value of Q when there is a sharp transition to a 
mixed phase. 

When we diminish the value of z 0: Fig.l, the behavior of eqc changes qualitatively. When 
Q/hu z is small the system behaves as in the previous case. However, when Q/hu z increases 



III. NUMERICAL RESULTS 
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there is a critical value of f2, f2 c<! where two branches of eqc appear, one stable, the other 
unstable. The two stable eqc distinguished by the separation of the centers of mass are 
called, respectively, distant and near stable eqc. When we further increase the value of Q, 
the relative distance of the centers of mass in the unstable eqc increases and merges with 
the stable distant eqc at a critical value of Q, Q c >, such that at Q > Q c> one is left with 
only the near stable eqc, where the condensates are mixed. For smaller values of z we have 
the same pattern, the values of Q c< where we have three branches of eqc and of Q c> where 
we have the merger of the unstable and distant stable eqc diminishing, this effect being less 
pronounced for the latter. 

Also shown in Fig.l is the graph of the branches of eqc for z = 0. We see that already 
at Q near zero we have three branches of eqc, where now in the near stable eqc there is 
complete overlap between the condensates (^> oa (z) = ip b(z)). When f2 increases, the totally 
mixed eqc remains, with a density profile independent of Q, whereas the distant stable and 
the unstable eqc approach each other and merge at f2 c> , such that, at Q > Q c> one is left 
with only the totally mixed eqc. Our results also show that we have spontaneous symmetry 
breaking effects at zq = [6,8]. Indeed, at zq = 0, our equations are separately invariant by 
space reflection, P z , and atom exchange P ex . The totally mixed eqc obey these symmetries 
separately, whereas the distant stable eqc do not, being invariant only by the product of 
these transformations. 

To find the signature of the onset of criticality, we calculate the normal modes along the 
branches of eqc. We can group the five normal modes into two sets. In one set there are two 
normal modes which are a linear combination of an out-of-phase oscillation of the centers 
of mass of each condensate and an in-phase breathing oscillation of the condensate densities 
with the center of mass of the mixture and the population fraction at its equilibrium values. 
In the second set we have three normal modes which are a linear combination of an in- 
phase oscillation of the centers of mass, out-of-phase breathing oscillation of the condensate 
densities and particle exchange between the condensates. The splitting of the normal modes 
into these two groups is a general result since it follows from the invariance of the equations 
(4-5) under P z P ex , the normal modes of the first group being even under this transformation 
and the one of the second, odd. 

We found that the signature of criticality in all cases is the collapse of a normal mode 
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which is an out-of-phase oscillation of center of mass (dipole oscillation) of the condensates 
and an in-phase breathing oscillation of the condensate densities. In fig. 2 we have a graph of 
the energies of the two normal modes of the first group, which involves the dipole oscillation 
of the condensates, for zq = 0.0Sz S h o , which is a case where we found the existence of critical 
points (see Fig.l). 

In Fig. 2a and Fig.2b we present a graph of the energies of the two normal modes along 
the distant and near stable eqc. In the distant eqc we see that there is one normal mode 
with an almost constant energy, Fig. 2a, and other whose energy increases beginning from 
fl = 0, reaches a maximum value and starts to decrease, Fig.2b. When fl approaches the 
critical value fl c >, where the distant stable eqc disappears, the energies of the two normal 
modes approaches each other and at fl = fl c> one of the energies goes very abruptly to zero. 
We have a similar behavior in the graph of the normal mode energies along the near stable 
eqc, Figs.2a-2b. Approaching from above the point where the near stable eqc disappears, 
the two energies approaches each other and again one the of them goes abruptly to zero at 
fl = fl c< . This behavior is completely general near a critical point. 

A "scar" of this critical behavior is also present for a value of zq at the interface between 
values where we have and we do not have critical points, such as zq = 0.23^0 (see Fig.l). 
We see in Figs.2c-2d that, corresponding to the very narrow range of values where the eqc 
change from separated to mixed, we have also an abrupt change in the values of the two 
normal modes energies. In Fig. 3a we detached the region of the sharp change and we see 
that it occurs in a very narrow range of values of fl and it is a consequence of a strong level 
repulsion between the two approaching even normal mode energies. 

In Fig. 3b we illustrate a generic phenomenon that occurs when fl — > 0, the appearance of 
a Goldstone zero energy mode. Indeed, when fl — > 0, the particle number of each component 
of the mixture is a conserved quantity and since our theory conserves only the total number 
of atoms, this violation is translated into the appearance of a zero energy mode. Fig. 3b 
shows how the energy of one of odd normal mode goes to zero for zq = 0.23z s h o - 

One question left untouched up to now is the identification of the lowest energy config- 
uration when we have many branches. In our model the answer is that, for small fl, the 
distant eqc is always the lowest energy configuration, changing to the near eqc at a higher 
value of fl, smaller than fl c> . However, the energy differences are very small, the equilibrium 
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configurations are almost degenerate. 

Our conclusions are based in calculations where we took \ aa = X bb = \ ab . It is well known 
that for homogeneous condensate mixtures, the parameters that controls the mechanism of 
phase separation is A aa Afef, — \ 2 ab [6]. In coupled mixtures, we add two additional factors which 
have opposite effects in the mechanism of spatial separation, the trap displacement and the 
laser coupling field and a point that deserves investigation is how robust are our conclusions 
when we relax the equal interaction strength condition. 

IV. SUMMARY 

To summarize, we have studied the equilibrium and stability properties of a coupled 
two-component BEC, as function of the laser field strength and trap displacement, using 
the variational method and the one-dimensional model of reference [14]. The laser field 
has a stabilizer effect in the mechanism of spatial separation of components in the mixture, 
opposite to the effect of the trap displacement. We found many branches of eqc, with a host 
of critical points. In all cases the signature of the onset of criticality is the collapse of a 
normal mode, which is a linear combination of an out-of-phase translation and an in-phase 
breathing oscillation of the condensate densities. 

When the traps are not displaced, we found eqc which exhibits symmetry breaking effects. 
In principle these eqc with broken symmetry can be reached by, starting at a sufficiently high 
value of Vt and z , adiabatically take the limit Q — > and z — > 0. Taking the limit in the 
opposite order, we end up in the symmetric eqc (see Fig.l). 

Undoubtedly our calculations are simple. However, it unveils a very rich structure in 
systems of coupled condensates, which should be explored experimentally and theoretically 
by more complete calculations. 
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Figure 1 
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Figure 1. The plot shows the relative distance of the centers of mass of the two condensates in 
the eqc, 2q la , as function of the laser strength f2, for fixed values of z . Curves from the top 
correspond to z = 0.3,0.23,0.08,0 z sho . The dotted and dashed curves indicate, respectively, 
the distant and near stable eqc. For z = the straight line 2g i(I = correspond to one branch 
of eqc. The laser strength is expressed in units of hw z and the distance in units of z s ho- See text 
for more details. 
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Figure 2 
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Figure 2. Fig. 2a and Fig. 2b show the energies of the two normal modes, which are an out-of-phase 
translation and an in-phase breathing oscillation of the condensate densities. The energies are 
calculated along the z = 0.08z sho curve, the dotted and the dashed lines correspond, respectively, 
to energies along the distant and near stable eqc. In Figs. 2c and 2d., we have a similar graph, 
now along the z = 0.23z sho curve. The energies are measured in units of Tiw z . See text for more 
details. 
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Figure 3 
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Figure 3. In Fig. 3a we have a plot of the energies near the point of sharp change, shown in 
Figs. 2c and 2d. In Fig. 3b we have a plot which shows how the energy of one of the odd normal 
modes goes to zero when f2 — > 0, for z = 0.23z zsho . See text for more details. 
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